home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / dgemv.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  7.4 KB  |  265 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
  5.      $                   BETA, Y, INCY )
  6. *     .. Scalar Arguments ..
  7.       DOUBLE PRECISION   ALPHA, BETA
  8.       INTEGER            INCX, INCY, LDA, M, N
  9.       CHARACTER*1        TRANS
  10. *     .. Array Arguments ..
  11.       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  DGEMV  performs one of the matrix-vector operations
  18. *
  19. *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
  20. *
  21. *  where alpha and beta are scalars, x and y are vectors and A is an
  22. *  m by n matrix.
  23. *
  24. *  Parameters
  25. *  ==========
  26. *
  27. *  TRANS  - CHARACTER*1.
  28. *           On entry, TRANS specifies the operation to be performed as
  29. *           follows:
  30. *
  31. *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
  32. *
  33. *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
  34. *
  35. *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  M      - INTEGER.
  40. *           On entry, M specifies the number of rows of the matrix A.
  41. *           M must be at least zero.
  42. *           Unchanged on exit.
  43. *
  44. *  N      - INTEGER.
  45. *           On entry, N specifies the number of columns of the matrix A.
  46. *           N must be at least zero.
  47. *           Unchanged on exit.
  48. *
  49. *  ALPHA  - DOUBLE PRECISION.
  50. *           On entry, ALPHA specifies the scalar alpha.
  51. *           Unchanged on exit.
  52. *
  53. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  54. *           Before entry, the leading m by n part of the array A must
  55. *           contain the matrix of coefficients.
  56. *           Unchanged on exit.
  57. *
  58. *  LDA    - INTEGER.
  59. *           On entry, LDA specifies the first dimension of A as declared
  60. *           in the calling (sub) program. LDA must be at least
  61. *           max( 1, m ).
  62. *           Unchanged on exit.
  63. *
  64. *  X      - DOUBLE PRECISION array of DIMENSION at least
  65. *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  66. *           and at least
  67. *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  68. *           Before entry, the incremented array X must contain the
  69. *           vector x.
  70. *           Unchanged on exit.
  71. *
  72. *  INCX   - INTEGER.
  73. *           On entry, INCX specifies the increment for the elements of
  74. *           X. INCX must not be zero.
  75. *           Unchanged on exit.
  76. *
  77. *  BETA   - DOUBLE PRECISION.
  78. *           On entry, BETA specifies the scalar beta. When BETA is
  79. *           supplied as zero then Y need not be set on input.
  80. *           Unchanged on exit.
  81. *
  82. *  Y      - DOUBLE PRECISION array of DIMENSION at least
  83. *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  84. *           and at least
  85. *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  86. *           Before entry with BETA non-zero, the incremented array Y
  87. *           must contain the vector y. On exit, Y is overwritten by the
  88. *           updated vector y.
  89. *
  90. *  INCY   - INTEGER.
  91. *           On entry, INCY specifies the increment for the elements of
  92. *           Y. INCY must not be zero.
  93. *           Unchanged on exit.
  94. *
  95. *
  96. *  Level 2 Blas routine.
  97. *
  98. *  -- Written on 22-October-1986.
  99. *     Jack Dongarra, Argonne National Lab.
  100. *     Jeremy Du Croz, Nag Central Office.
  101. *     Sven Hammarling, Nag Central Office.
  102. *     Richard Hanson, Sandia National Labs.
  103. *
  104. *
  105. *     .. Parameters ..
  106.       DOUBLE PRECISION   ONE         , ZERO
  107.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  108. *     .. Local Scalars ..
  109.       DOUBLE PRECISION   TEMP
  110.       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
  111. *     .. External Functions ..
  112.       LOGICAL            LSAME
  113.       EXTERNAL           LSAME
  114. *     .. External Subroutines ..
  115.       EXTERNAL           XERBLA
  116. *     .. Intrinsic Functions ..
  117.       INTRINSIC          MAX
  118. *     ..
  119. *     .. Executable Statements ..
  120. *
  121. *     Test the input parameters.
  122. *
  123.       INFO = 0
  124.       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
  125.      $         .NOT.LSAME( TRANS, 'T' ).AND.
  126.      $         .NOT.LSAME( TRANS, 'C' )      )THEN
  127.          INFO = 1
  128.       ELSE IF( M.LT.0 )THEN
  129.          INFO = 2
  130.       ELSE IF( N.LT.0 )THEN
  131.          INFO = 3
  132.       ELSE IF( LDA.LT.MAX( 1, M ) )THEN
  133.          INFO = 6
  134.       ELSE IF( INCX.EQ.0 )THEN
  135.          INFO = 8
  136.       ELSE IF( INCY.EQ.0 )THEN
  137.          INFO = 11
  138.       END IF
  139.       IF( INFO.NE.0 )THEN
  140.          CALL XERBLA( 'DGEMV ', INFO )
  141.          RETURN
  142.       END IF
  143. *
  144. *     Quick return if possible.
  145. *
  146.       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  147.      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  148.      $   RETURN
  149. *
  150. *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
  151. *     up the start points in  X  and  Y.
  152. *
  153.       IF( LSAME( TRANS, 'N' ) )THEN
  154.          LENX = N
  155.          LENY = M
  156.       ELSE
  157.          LENX = M
  158.          LENY = N
  159.       END IF
  160.       IF( INCX.GT.0 )THEN
  161.          KX = 1
  162.       ELSE
  163.          KX = 1 - ( LENX - 1 )*INCX
  164.       END IF
  165.       IF( INCY.GT.0 )THEN
  166.          KY = 1
  167.       ELSE
  168.          KY = 1 - ( LENY - 1 )*INCY
  169.       END IF
  170. *
  171. *     Start the operations. In this version the elements of A are
  172. *     accessed sequentially with one pass through A.
  173. *
  174. *     First form  y := beta*y.
  175. *
  176.       IF( BETA.NE.ONE )THEN
  177.          IF( INCY.EQ.1 )THEN
  178.             IF( BETA.EQ.ZERO )THEN
  179.                DO 10, I = 1, LENY
  180.                   Y( I ) = ZERO
  181.    10          CONTINUE
  182.             ELSE
  183.                DO 20, I = 1, LENY
  184.                   Y( I ) = BETA*Y( I )
  185.    20          CONTINUE
  186.             END IF
  187.          ELSE
  188.             IY = KY
  189.             IF( BETA.EQ.ZERO )THEN
  190.                DO 30, I = 1, LENY
  191.                   Y( IY ) = ZERO
  192.                   IY      = IY   + INCY
  193.    30          CONTINUE
  194.             ELSE
  195.                DO 40, I = 1, LENY
  196.                   Y( IY ) = BETA*Y( IY )
  197.                   IY      = IY           + INCY
  198.    40          CONTINUE
  199.             END IF
  200.          END IF
  201.       END IF
  202.       IF( ALPHA.EQ.ZERO )
  203.      $   RETURN
  204.       IF( LSAME( TRANS, 'N' ) )THEN
  205. *
  206. *        Form  y := alpha*A*x + y.
  207. *
  208.          JX = KX
  209.          IF( INCY.EQ.1 )THEN
  210.             DO 60, J = 1, N
  211.                IF( X( JX ).NE.ZERO )THEN
  212.                   TEMP = ALPHA*X( JX )
  213.                   DO 50, I = 1, M
  214.                      Y( I ) = Y( I ) + TEMP*A( I, J )
  215.    50             CONTINUE
  216.                END IF
  217.                JX = JX + INCX
  218.    60       CONTINUE
  219.          ELSE
  220.             DO 80, J = 1, N
  221.                IF( X( JX ).NE.ZERO )THEN
  222.                   TEMP = ALPHA*X( JX )
  223.                   IY   = KY
  224.                   DO 70, I = 1, M
  225.                      Y( IY ) = Y( IY ) + TEMP*A( I, J )
  226.                      IY      = IY      + INCY
  227.    70             CONTINUE
  228.                END IF
  229.                JX = JX + INCX
  230.    80       CONTINUE
  231.          END IF
  232.       ELSE
  233. *
  234. *        Form  y := alpha*A'*x + y.
  235. *
  236.          JY = KY
  237.          IF( INCX.EQ.1 )THEN
  238.             DO 100, J = 1, N
  239.                TEMP = ZERO
  240.                DO 90, I = 1, M
  241.                   TEMP = TEMP + A( I, J )*X( I )
  242.    90          CONTINUE
  243.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  244.                JY      = JY      + INCY
  245.   100       CONTINUE
  246.          ELSE
  247.             DO 120, J = 1, N
  248.                TEMP = ZERO
  249.                IX   = KX
  250.                DO 110, I = 1, M
  251.                   TEMP = TEMP + A( I, J )*X( IX )
  252.                   IX   = IX   + INCX
  253.   110          CONTINUE
  254.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  255.                JY      = JY      + INCY
  256.   120       CONTINUE
  257.          END IF
  258.       END IF
  259. *
  260.       RETURN
  261. *
  262. *     End of DGEMV .
  263. *
  264.       END
  265.